Dynamics of channel incision in a granular bed driven by subsurface water flow 
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We propose a dynamical model for the erosive growth of a channel in a granular medium driven 
by subsurface water flow. The model is inferred from experimental data acquired with a laser-aided 
imaging technique. The evolution equation for transverse sections of a channel has the form of a 
non-locally driven Burgers equation. With fixed coefficients this equation admits an asymptotic 
similarity solution. Ratios of the granular transport coefficients can therefore be extracted from the 
shape of channels that have evolved in steady driving conditions. 
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One of the salient features of fluvial erosion is the 
formation and growth of channels 1]. The channels 
act to focus the flow of water, either over land 0, 
or, ubiquitously but less commonly appreciated, be- 
neath the surface Q. Here we address the dynam- 
ics of erosion driven by subsurface flow. Whereas flow 
through the subsurface is well characterized by Darcy's 
law 0, , channel formation and growth resulting from 
seepage flows out of the surface is relatively poorly un- 
derstood 0. This erosive process is important not only 
because of its relevance to enigmatic features of natural 
landscapes 0, S SH, E3, EH EE E3, E3 , but also because 
it raises fundamental questions concer ning the continuum 
mechanics of wet sand [llllTrilHITllirlEolEllE^I. 

In this Letter we construct a theoretical model of the 
erosive dynamics of an isolated channel. The model is 
inferred from data acquired with a laser-aided imaging 
technique. In our experiment we grow a single channel 
under controlled driving conditions and acquire maps of 
the evolving granular surface that are highly resolved in 
space and time. We find that transverse sections of a 
channel evolve according to a non-locally driven Burgers 
equation. Under steady driving conditions the evolution 
equation admits a similarity solution, thereby relating 
erosive transport coefficients to channel geometry. 
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FIG. 1: (a) Schematic of our experimental setup, (b) Example 
of an isolated channel observed after 90 minutes of subsurface 
water flow driven growth. A small initial channel is incised 
on the surface to fix the location and number of channels. 
Similar shapes are observed if channels form spontaneously. 

Figure shows a schematic of the experiment, as- 
pects of which are described in more detail elsewhere 



[fj, |2j|. Water enters beneath a pile of identical glass 
beads through a fine mesh and exits at the foot of the 
pile through the same kind of mesh. The water flux is 
controlled by the height H of the water column in a reser- 
voir behind the pile. The slope of the initial sandpilc 
as well as the water column height are the control vari- 
ables of the experiment. Since the water table is con- 
vex upward in this geometry |23| . we are able to finely 
control the amount of water on the surface of the gran- 
ular pile. Refs. and describe the phenomenology 
of the pattern formation in this setup and quantify the 
transitions between various modes of granular flow: sur- 
face flow (responsible for the formation of the channel 
network), slumping (bulk frictional instability) and flu- 
idization. 

Here we focus on the late-stage evolution of an isolated 
channel that grows from a short initial channel of trian- 
gular crossection. A scanning laser imaging technique 
allows us to measure the evolving height of the sandpilc 
with sub-grain resolution in space and one-minute res- 
olution in time. A laser sheet scans the surface while 
a digital camera acquires images from an oblique angle. 
The height of the surface is then extracted from an im- 
age of the intersection of the laser sheet with the granular 
surface. 

The water level H is set below the threshold for erosion 
outside the channel but above the threshold for erosion 
in the channel. An example of a channel grown from an 
initial channel with a triangular crossection is shown in 
Fig. IHb) . We find that the late stage morphology of the 
channel is insensitive to the exact initial condition as long 
as the initial incision is sufficiently deep. 

We seek an effective description of the channel's evo- 
lution in terms of the surface height h(x, y, t) measured 
from the uneroded surface. In other words, we seek an 
expression for the erosion rate h = dh/dt in terms of 
h, its spatial derivatives, and possibly spatial integrals. 
(Here y is the downslope axis and x is the axis transverse 
to the channel). This approach is reasonable because the 
relaxation of the water flux toward steady state is fast 
compared to the erosion rate. Thus the water flux is a 
functional of the slowly changing shape of the channel. 
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FIG. 2: Residual erosion rate at t = 39 min plotted against (a) the slope, (b) the curvature, and (c) the fractional channel 
depth. Each symbol represents an average over several hundred data points. To obtain the residual erosion rate we subtract 
from the measured h the terms on the right hand side of Eq. Q which depend on the variables other than the one against 
which the erosion rate is being plotted. Lines are drawn using the parameters extracted by a least squares fit of the unbinned 
data to equation Q. Error bars represent two standard deviations of the mean. 



The erosion rate is a function of the local water and 
granular fluxes. In steady state these fluxes are function- 
als of the global shape of the granular pile. We seek to 
reduce the global dependence to a single scalar by argu- 
ing that, because the channel evolves in a roughly self- 
similar manner, the water fluxes (bulk as well as surficial) 
are functions of the local topography and a time- varying 
scale factor. We choose the channel depth h (y) as this 
scale factor (h corresponds to height of the deepest point 
of the crossection) . We have thus reduced the problem 
to finding the dependence of the erosion rate on the local 
topography and the overall global scale factor. 

We make one more simplification: we consider solely 
the evolution of transverse sections of the channel, and 
therefore express the erosion rate h = dh/dt through h 
and its derivatives with respect to the transverse coordi- 
nate x only. Variation of the water flow in the downs- 
lope y-direction is accounted for via y-dependent coef- 
ficients. This approximation is reasonable everywhere 
except the top portion of the channel's head where the 
downslope gradient varies rapidly. Grains and micro- 
scopic avalanches enter and leave a given channel tran- 
sect. Projected onto this transect, the transport of height 
h is no longer volume conserving. 

Thus, the erosion rate h is assumed to be a function 
of the fractional depth h/h®, the transverse slope h x — 
dh/dx, and the curvature h xx — d 2 h/dx 2 . Our approach 
is to measure these quantities for data in a window of 
duration At in time and length Ay in the downslope 
coordinate and to fit the resulting data cloud via a least 
squares method to the form 

h = vh xx - b\h x \ -Xh 2 x -fi Q(h/h - /), (1) 

where O(x) is the Heaviside step function. The empirical 
constants /i, v, A, 5, and / are functions of time t and 
downslope coordinate y. They encode the microscopic 
properties of the grain dynamics as well as the strength 
of the driving water flow. The diffusion constant v re- 



flects the rate of smoothing of local perturbations. The 
last term on the right hand side represents driving due 
to the seeping water. We assume that driving is constant 
below a fraction / of the total channel depth and null 
above this fractional depth. The second and third terms 
on the right hand side of are advective. We hypoth- 
esize that perturbations are advected only up the slope, 
similar to Ref. 24] . Therefore, the non-analytic term 
S\h x \ corresponds to advection of perturbations with ve- 
locity 5 independent of slope, whereas Xh 2 corresponds 
to advection with velocity Xh x , which grows linearly with 
slope. 

Figurc|2illustrates the fit of the data to the model. As 
shown in Fig. Gla), the slope dependent terms make the 
largest contribution into the erosion rate. Data points 
with slopes smaller than 0.7, comprising over 95% of all 
data, have been used in the fit. A discrepancy between 
the model and the data occurs for larger slopes where the 
erosion rate saturates. A measurable contribution to the 
erosion rate due to diffusion is shown in Fig. Elb) . The 
extracted positive granular diffusivity ;/ is statistically 
significant. The driving, depth-dependent term, shown 
in Fig.[5Jc), is approximated by a step function although 
the data show a more gradual transition. We expect the 
width of the transition region to remain constant as the 
channel grows. Thus the step function approximation 
should improve with time. 

We extract the timc-dcpcndcnt coefficients in Eq. 
near a transect fixed in the lab frame. Fig. shows the 
extracted parameters for a transect initially above the 
pre-dug channel. As the channel head migrates past 
this transect, the driving fj, peaks and declines. The 
values of advection speeds 6 and A also decrease since 
they are related to the driving. The diffusion coefficient 
v = 0.039 ± 0.033 mm 2 /min and the fractional driving 
/ = 0.52 ± 0.14 fluctuate around their respective means 
which are roughly time- independent. 

In Fig.0|a) and (b) we indicate in a contour plot of the 
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FIG. 3: Example of the evolution of the coefficients (i, A and S 
(measured in mm/min) extracted using Eq. Q from the data 
using taken near a fixed transect illustrated in Fig. 2] The 
head of the channel arrives to this transect approximately one 
minute after the start of the experimental run. 



channel the location of the transect for which the param- 
eters in Fig.|21are computed. Fig.QJc) compares the mea- 
sured shapes of this transect with the shapes evolved via 
Eq. with the time-dependent coefficients. Agreement 
indicates that the right hand side of Eq. JIJ captures the 
essential features of the erosion rate in the channel ge- 
ometry. The rate of advance of the channel sidewalls and 
the receding rate of channel bottom are determined by 
different combinations of the parameters. The match of 
the evolved shapes to the data shows that both rates are 
in accord with the experiment. 

In our experiment the coefficients, intimately related to 
the geometry of the water flow, change markedly during 
the evolution of the channel as illustrated in Fig. [3J It is 
conceivable, however, that in other geometries the trans- 
port coefficients are approximately constant over a long 
time. It is therefore appropriate to examine the long-time 
behavior of Eq. with constant coefficients. The rest 
of the paper is devoted to the asymptotic calculation. 

The main result is that the non-local driving term in 
Eq. QJ allows for an asymptotically self-similar growing 
channel shape. To see that let us define a scale fac- 
tor £(t) and the shape Y(Z> t), which is a function of 
the scaled coordinate £ = xjl(t) via a transformation 
h(x, t) = £(t) Y(Z, t). Further progress results from a for- 
mal expansion 

Y(U) = Y (O + ]y 1 (O + .-., (2a) 
t{t) = t + £ 1 lnt + .... (2b) 

The zeroth order shape Yq(Z) is independent of time. 
Thus, if the expansion J2J converges, the shape converges 
(albeit slowly as 1/t) to a similarity solution Yq. A sub- 
stitution of the expansion J5J into QJ yields, at the lowest 
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FIG. 4: Contour plots of the laser height data at (a) t = 
15 min and (b) 115 min after the start of the water flow. 
The thick line shows the location of the modeled transect. 
The contour interval is 5 mm. (c) Same transect at t = 15, 
45, and 115 min compared to the prediction of Equation Jy. 
The smoothed shape of the transect at t = 15 min was used 
as the initial condition. Actual time dependent parameters 
presented in Fig. [3]were used. 

order, 

Yo - £F ' = -PYq - a(*o) 2 - © (yr|j - /) » (3) 

where primes denote differentiation with respect to £. We 
scaled lengths by a — v/n and time by t = vj ^ and de- 
fined = and a = A/ fi. Note that the diffusive term 
does not enter at the lowest order. The physical reason 
for this is that, as we show below, diffusion is important 
on a fixed length scale. As the channel grows larger than 
this diffusive length, the h xx term in JIJ becomes negli- 
gible compared to the advection and the driving terms. 

The solution to J3J must be symmetric with respect 
to £ — * — £ and smooth everywhere except at the frac- 
tional driving point Zo where Yo(£o) = /Y)(0). It turns 
out that there exists a one-parameter family of similarity 
solutions which satisfy these criteria. These solutions are 
constructed as follows. In the driven region any piece of 
the parabola 

Yo(0 = ^(£ 2 -2/^ + /? 2 -4a), (4) 

and a tangent line to this parabola at a point £i, 

r (0 = ^(2£(£i"/3)+/3 2 -£ 1 2 -4a), (5) 

are solutions to J3J- In the undriven region, the trivial 
solution ft = as well as tangents to the parabola 

Y (o = ^(e-m) (6) 

are solutions. A smooth solution symmetric around £ = 
is therefore constructed from four pieces: 

0<£<£ b , 

+ ei <£<£ 

-f + Soit-to), Zo<Z<Ze- 
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The first piece is the flat bottom of the channel tangent 
to the parabola Q at its apex = 0, Yo = — 1. The 
second piece is part of the parabola itself. The third piece 
is another tangent to parabola Q at a point £i such that 
P < £i < /3+\/4a(l — /)■ The fourth piece is the tangent 
of slope so to the undriven parabola ©■ Note that Yq is 
not smooth at £ (driving) and £ e (channel edge). 

It turns out that the diffusive term, though not present 
at zeroth order, acts to select a unique member from the 
one-parameter family of similarity solutions. The selec- 
tion mechanism, verified by numerical methods, is un- 
clear to us at this time. The selected similarity solution 
is one in which the second tangent to the parabola |@J 
is missing, i.e., £i assumes its upper limit. The slope of 
the sidewall of the channel in the undriven region is then 
simply 



SQ = —(l + ^]—f). (8) 

The expressions for £o and £ e are simple as well: 

Co = P + v / MW), £e = £o + f/s (9) 

We remark that the asymptotic shape consisting of a 
flat bottom, curved parabolic flanks in the driven region 
followed by straight sidewalls can be used successfully to 
fit channel crossections which evolve in non-steady con- 
ditions. However, the ratios of parameters extracted by 
such an asymptotic fit can deviate greatly from the pa- 
rameters extracted by the fit of the equation to the 
data cloud, because the shape of the channel at the time 
of the fit retains a memory of its prior dynamical state. 

Besides selecting the unique self-similar shape, the dif- 
fusive term also acts to smooth slope discontinuities at 
£o and £ e . An exact expression for this smoothing can be 
obtained at the channel's edge £ e . In the vicinity of this 
point, the slope s = h x is a hyperbolic tangent 

S (M) = y(l-tanh^ (£-«*)), (10) 

moving with velocity v — sqcy + (3 and smooth on a scale 



As we mentioned above, diffusion acts on a fixed scale R. 
The smoothing of the kink at £ = £o is probably related 
to the selection of the unique similarity solution. 

Guided by an experiment, we have developed a phe- 
nomenological model of granular dynamics in a transect 
of a channel eroded by subsurface fluid flow. Precise, 
time resolved measurements of the height of the eroding 
sandpile allow us to test the validity of the model. The 
erosion rate is a function of the local topography and an 
overall scale factor. We expect this approximation to be 



good as long as the channel's shape remains roughly self- 
similar. There are five parameters in the model which in 
our experiment vary with time. These parameters encode 
the geometry of the water flow as well as the features of 
the microscopic granular flow. In a different water flow 
geometry these parameters can remain roughly constant. 
In this case the height evolution equation admits a 
self-similar solution and the three dimensionlcss parame- 
ters a, j3 and / dictate the self-similar shape. Conversely, 
if a shape is known to have resulted from evolution with 
roughly constant parameters, then a, (3 and / can be 
extracted from the static shape. Given an independent 
measure of the erosion rate [i at the bottom of the chan- 
nel, as well as the measurement of the smoothing length 
scale R at the channel's edge, all five parameters can 
be recovered. Such a procedure should be useful for the 
analysis of natural channels formed by subsurface fluid 
flow. 

This work is supported by DOE grants DE- 
FG02-02ER15367 at Clark University and DE-FG02- 
99ER15004 at MIT. 



[1] W. E. Dietrich and T. Dunne, in Channel Network Hy- 
drology, edited by K. Beven and M. J. Kirby (John Wiley 
and Sons Ltd, New York, 1993), pp. 175-219. 

[2] R. E. Horton, Geological Society of America Bulletin 56, 
275 (1945). 

[3] T. Dunne, in Special Paper 252 (Geological Society of 
America, Boulder, CO, 1990), pp. 1-28. 

[4] A. E. Sheidegger, The Physics of Flow Through Porous 
Media (Macmillan Company, New York, 1960). 

[5] J. Bear, Dynamics of Fluids in Porous Media (Dover 
Publications, New York, 1972). 

[6] N. Schorghofer, B. Jensen, A. Kudrolli, and D. H. Roth- 
man, J. Fluid Mech. 503, 357 (2004). 

[7] A. D. Howard and C. F. McLane, Water Resources Re- 
search 24, 1659 (1988). 

[8] J. E. Laity and M. C. Malin, Geol. Soc. Am. Bull. 96, 
203 (1985). 

[9] S. A. Schumm, K. F. Boyd, C. G. Wolff, and W. J. Spitz, 

Geomorphology 12, 281 (1995). 
[10] D. L. Orange, R. S. Anderson, and N. A. Breen, GSA 

Today 4, 1 (1994). 
[11] C. G. Higgins, Geology 10, 147 (1982). 
[12] C. K. Wentworth, J. Geol. 36, 385 (1928). 
[13] R. C. Kochel and J. F. Piper, J. Geophys. Res. 91, E175 

(1986). 

[14] M. C. Malin and K. Edgett, Science 288, 2330 (2000). 
[15] N. Huang, G. Ovarlez, F. Bertrand, S. Rodts, P. Coussot, 

and D. Bonn, Phys. Rev. Lett. 94, 028301 (2005). 
[16] Z. Fournier, D. Geromichalos, S. Herminghaus, M. M. 

Kohonen, F. Mugele, M. Scheel, M. Schulz, B. Schulz, 

C. Schier, R. Seemann, et al., J. Phys.-Cond. Mat. 17, 

S477 (2005). 

[17] P. Tegzes, T. Vicsek, and P. Schiffer, Phys. Rev. E 67, 
051303 (2003). 

[18] M. Schulz, B. M. Schulz, and S. Herminghaus, Phys. Rev. 



■5 



E 67, 052301 (2003). 

[19] N. Jain, J. M. Ottino, and R. M. Lueptow, J. FLuid 
Mech. 508, 23 (2004). 

[20] N. Jain, D. V. Khakhar, R. M. Lueptow, and J. M. Ot- 
tino, Phys. Rev. E 86, 3771 (2001). 

[21] J. Zhou, A. Bertozzi, B. Dupuy, and A. E. Hosoi, Phys. 
Rev. Lett. 94, 117803 (2005). 



[22] A. Daerr, P. Lee, J. Lanuza, and E. Clement, Phys. Rev. 

E 67, 065201 (pages 4) (2003). 
[23] A. E. Lobkovsky, W. Jensen, A. Kudrolli, and D. H. 

Rothman, J. Geophys. Res-Earth Surface 109, 1 (2004). 
[24] T. Boutreux, E. Raphael, and P. G. de Gennes, Phys 

Rev. E 58, 4692 (1998). 



